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A microfluidic device is constructed from PDMS with a single channel having a 

short section that is a thin flexible membrane, in order to investigate the complex 

fluid-structure interaction that arises between a flowing fluid and a deformable 

wall. Experimental measurements of membrane deformation and pressure drop 

are compared with predictions of two-dimensional and three-dimensional compu- 
te ', 

tational models which numerically solve the equations governing the elasticity of 

the membrane coupled with the equations of motion for the fluid. It is shown that 

the two-dimensional model, which assumes a finite thickness elastic beam that is 

infinitely wide, approximates reasonably well the three-dimensional model, and is 

in excellent agreement with experimental observations of the profile of the mem- 

brane, when the width of the membrane is beyond a critical thickness, determined 

\Q ', to be roughly twice the length of the membrane. 
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I. INTRODUCTION 

The fabrication of microfluidic devices from soft polymers or elastomers has gained 
considerable interest in the last decade. The attractiveness of using these soft materials, 
in particular, stems from the ability to tailor the polymer's physicochemical properties 
specifically for a given application, the lower material and fabrication costs which allows 
the possibility for disposable devices, and the durability of polymer-based materials 
compared to the brittleness of conventional hard materials such as silicon and glass 1 . 
Moreover, soft polymers such as polydimethysiloxane (PDMS) 2 offer excellent optical 
transparency, gas permeability and biocompatibility, vital for on-chip cell culture and 
comprising a large proportion of microfluidic devices for cellomics, drug screening and 
tissue engineering 3 . 

The bonding strength, ability to mold at even nano-scale, biocompatibility, trans- 
parency and flexibility of these elastomeric substrates also make them ideal for fabricat- 
ing microfluidic actuation structures. For example, thin PDMS membranes have been 
employed as diaphragms or membrane interfaces for pneumatic actuation and control 
in microchannels 4 " 7 and substrates for biological characterisation and manipulation in 
microdevices 8 " 10 More sophisticated microfluidic actuation structures have also been pro- 
posed, including multilayer and branched channel networks controlled by elastomeric 
micropumps and microvalves 11 for a variety of uses, for example, to spatiotemporally 
control chemical gradients for chemotaxis studies on a microfluidic chip 12 . 

Fundamental studies in investigating the complex fluid-structural interaction arising 
from these flexible materials and the flow of the fluid within them, however, have not 
kept at the same rapid pace as the developments that have arisen. In fact, there have 
been no studies undertaken to investigate the flow through flexible channels at scales 
commensurate with microfluidic devices. To the best of our knowledge, even simple 
experiments 13 " 22 and theoretical studies 23 " 28 undertaken to investigate Newtonian flows 
through macroscopic deformable tubes have yet to be reproduced at the microscale, where 
channel dimensions are of the order of 10-100 [im and the Reynolds number is typically 
of the order of unity or below, typically two or more orders of magnitude smaller than 
the O(10 mm) channel dimensions and O(100) or greater Reynolds numbers examined in 
these studies. The length scale at which fluid flow occurs in microfluidic devices is entirely 
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different from the large-scale flows. Fluid flowing in a conventional microfluidic channel 
with characteristic length scale in the sub-millimeter range, is identified by low velocity 
and hence small Reynolds numbers. It is widely acknowledged that the experimental 
observations conducted in macroscale channels can be well predicted by the Navier-Stokes 
equation. However, to the best of our knowledge, there is no evidence in the literature of 
the use of a fluid-structural interaction theory to a collapsible microchannel. Fluid flow 
within a flexible structure is regulated by the stresses imposed upon the structure by both 
the fluid and any external forces. Thus, the rheological properties of both the fluid and 
the structure significantly influence the fluid flow in the system. The viscous stresses and 
fluid pressure exerted on the boundaries of the flexible wall cause its deformation. Due to 
the deformation of the flexible structures, the flow domain and flow field alters and gives 
rise to an intricate fluid-structure interaction problem which requires the solution of a free- 
boundary problem. Thus in order to understand fluid-structural interaction phenomenon 
at the microscale, successful development of fabricated microfluidic devices as well as the 
implementation of a fluid-structural interaction theory associated with the microscale is 
necessary. 

Given this motivation, we therefore carry out a fundamental investigation of the fluid 
flow in a deformable raz'crochannel. Specifically, we compare deformation profiles from 
experiments carried out on a custom fabricated microfludics with a flexible membrane 
section with that predicted by a two-dimensional finite element model that solves the 
coupling between the membrane deformation and the fluid flow. We fabricated the 
structure by casting PDMS into a 200 fim high and 29 mm long microchannel. Membrane 
deformation was controlled in the experiment by the level of air pressure introduced via 
a pressure chamber from the side of the membrane opposing the channel. 

The rest of the article is thus organised as follows. We first formulate the numerical 
model and discuss the solution methodology in Sec. II. The fabrication and design of the 
deformable microchannel and the experimental methodology is subsequently described 
in Sec. III. A comparison between the results obtained from both the experiments and 
numerical simulation then follows in Sec. IV, after which we summarise our conclusions 
in Sec. V. 
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FIG. 1. Schematic illustration of the experimental setup. 



II. NUMERICAL SIMULATION AND SOLUTION METHODOLOGY 

A. Two-Dimensional Finite Element Model for Fluid-structure Interaction 
(2D-FEM-FSI) 



To match the geometry of the microchannel and the flexible membrane that spans the 
channel width in the experimental design (discussed subsequently in Sec. Ill) shown in Fig. 
1, we consider a two-dimensional model of the experimental setup as shown in Fig. 2(a), 
in which fluid flows through a section of the microchannel with height H, along that, a 
short segment of elastic membrane BC with thickness t and length L spanning the width 
of the channel exists on one side. Whilst the sidewalls of the channel prior to and after the 
membrane section, AB and CD with lengths L u and L^, respectively, are considered rigid, 
the membrane is allowed to deform under an external pressure p e as measured in the air 
pressure chamber. To mimic the experimental geometry, we set H = 200 /.mi, L u = 14 mm, 
L = 1 mm, La = 14 mm and t = 60 um. Here, the x-axis spans the channel length whereas 
the z-axis denotes the height of the channel with origin at point O. 

In the absence of body forces, assumed negligible at the microscopic scales considered, 
the equations of motion for steady, incompressible flow are specified by the continuity 
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FIG. 2. (a) Geometry of the two-dimensional deformable microchannel defining the solution 
space for the numerical simulation, (b) The solution strategy for the coupling between the fluid 
and solid domains that involves a free boundary problem is handled by mapping the physical 
fluid domain to a reference computational domain (Qf — > Qqp) and the physical solid domain to 
a reference zero-stress domain (Qg — > Os)- We note that the solution of the elasticity equations 
itself constitutes a mapping from the zero-stress configuration (Qg — > £1$), and is consequently not 
solved separately. The mapping from computational domain (Cos) to the zero-stress configuration 
(Qg) is known and it only requires a change of domain of integration. 



and momentum conservation, respectively: 



V • v = 0, 

5 



(1) 
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p v . Vu = V • (-/? I + t), (2) 

where p is the density of the fluid, v is the liquid velocity field, p the liquid pressure and r 
the viscous stress tensor; I represents the identity tensor. For a Newtonian fluid, r = 2r/D, 
where r] is the viscosity of the liquid and D = \{Vv + Vv T ) is the strain rate tensor. 

Given the deformability of the membrane under external pressure, the system com- 
prises a free boundary problem in which the fluid flow and the solid domain that con- 
stitutes the membrane are coupled. As we are interested in the steady flow, the inertia 
of the solid component is not affecting the overall dynamics of the system. We therefore 
employ a solution strategy that continuously maps the fluid and solid domains x = x(£), 
both unknown a priori, to arbitrary reference domains following the approach of Carvalho 
and Scriven 29 , as illustrated schematically in Fig. 2(b). Here, the physical and reference 
computational domains are parameterised by the position vector x and £, respectively, 
and X represents the position in the reference stress-free domain. The physical fluid 
domain Of is mapped by elliptic mesh generation to a reference computational domain 
Qof, where Eqs. (1) and (2) are solved. Due to the complexity in the geometry, the physical 
domain cannot be mapped to a simpler, quadrangular reference domain. Instead, it is 
more convenient to subdivide the physical domain into subdomains and then map each 
subdomain into a separate subdomain of the reference computational domain. Here we 
use a boundary-fitted finite element based elliptic mesh generation method 30 " 33 which 
involves solving the following elliptic differential equation for the mapping: 

V • D • V£ = 0, (3) 

where the dyadic D is a function of £ in a manner analogous to a diffusion coefficient, 
which controls the spacing of the coordinate lines 32 . 

The physically deformed solid domain constituting the flexible membrane Qs, on 
the other hand, is mapped to a reference domain, that for convenience, we choose as a 
hypothetical zero-stress state where the stress tensor vanishes over the entire membrane 
(which may not and need not be physically realised). It is in this stress-free domain 
Q s where the elasticity equations governing the deformation of the solid membrane are 
solved, although the solution of these equations itself constitutes a mapping from the zero- 
stress configuration Q s to the deformed domain Q s . The mapping from the computational 
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domain (Qos) to the zero-stress configuration (Qg) is known and only requires a change 
of the domain of integration. 

In the reference stress-free domain Q s , the equilibrium equation that governs the solid 
if acceleration and body forces can be neglected, reads 

V x • S = 0, (4) 

where S is the first Piola-Kirchhoff stress tensor. We note that this is related to the original 
deformed state of the solid, i.e., the physical solid domain, through Piola's transformation 
to the Cauchy stress tensor a- by 

S = F 1 • <r, (5) 

where, 

is the deformation gradient tensor, which relates the undeformed state X = (X, Y, Z) to 
the deformed state x = (x, y, z). Closure to the above is obtained through a constitutive 
relationship that relates the Cauchy stress tensor with the strain. For a neo-Hookean 
material, this takes the form 

ct = -n* I + GB, (7) 

where n* is a pressure-like scalar function, G is the shear modulus and B = F • F T . B is the 
left Cauchy-Green deformation tensor. 

The equations governing the fluid motion and the solid deformation above are subject 
to the following boundary conditions: 

1. No slip boundary conditions apply on the rigid walls, i.e., v = on -(L u + L/2) < 
x<(L d + L/2) when z = and -(L u + L/2) <x< -L/2 and L/2 < x < (L d + L/2) when 
z = H. 

2. Zero displacements are prescribed at the left side and right side of the solid, i.e., x 
= X on H<z<H + t when x = -L/2, L/2. 

3. At the upstream boundary (0 < z < H, x = -(L u + L/2)), a fully-developed velocity 
profile is specified, i.e., v z = and v x = U f(z/H) where U is the average inlet 
velocity. 
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4. At the downstream boundary (0 < z < H, x = (L^ + L/2)), the fully-developed flow 
boundary condition is imposed, i.e., n • V« = 0. 

5. A force balance and a no-penetration condition are prescribed at the interface be- 
tween the liquid and solid domain: 

n • t = n • <j and u = v = 0, (8) 

where u is the velocity of the solid and n is the outward unit vector normal to the 
deformed solid surface. 

6. A force balance is prescribed at the top surface. 

n.cr = -p e n, (9) 

where p e is the external pressure. 

7. We have compared our theory with the pressure drop only. So, the gauge pressure 
of the fluid at the downstream boundary is assumed negligible, i.e.., pd = 0. 

The weighted residual form of Eqs. (l)-(4), obtained by multiplying the governing 
equations with appropriate weighting functions and subsequently integrating over the 
current domain, yields a large set of coupled non-linear algebraic equations, which is 
solved subject to the boundary conditions specified above using Newton's method with 
analytical Jacobian, frontal solver and first order arc length continuation in parame- 
ters 33 " 36 . The formulation of the fluid-structure interaction problem posed here follows 
the procedure introduced previously by Carvalho and Scriven 29 in their examination of 
roll cover deformation in coating flows. It turns out, however, that the weighted residual 
form of Eq. (4) used in their finite element formulation is incorrect. While insignificant 
for small deformations, this error leads to significant discrepancies when the deformation 
is large. The weighted-residual equation is corrected here and validated in Appendix A 
against predictions using a commercial finite element software ANSYS 11.0 37 for the 
deformation of a simple beam fixed at its edges that we describe next. 

8 
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B. Finite Element Model in ANSYS (2D/3D-ANSYS) 

We expect that the two-dimensional numerical simulation proposed above reasonably 
approximates a three-dimensional system when the microchannel and hence membrane 
width is large compared to the height and length of the microchannel such that boundary 
effects at the membrane edges can be neglected. To determine the limits of the microchan- 
nel width at which the two-dimensional approximation breaks down, we carry out a 
three-dimensional finite element simulation that involves a plane-strain model for a com- 
pressible neo-Hookean solid since the incompressible neo-Hookean model described in 
Sec. II A is only valid for a two-dimensional geometry. For a given strain-energy density 
function or a elastic potential function of a neo-Hookean material, W, 

W=|(I 1 -3) + i(/-l) 2 , (10) 

where I\ = tr(C) is the first invariant of the right Cauchy-Green deformation tensor, G 
the initial shear modulus of the material, d the material incompressibility parameter and 
/ = det(F) the ratio of the deformed elastic volume over the undeformed volume of 
material. The corresponding stress component is 

s = 2 lc' (11) 

where S is the second Piola-Kirchoff stress tensor and C = F T • F is the right Cauchy-Green 
strain tensor. If acceleration and body forces are negligible, the equilibrium equation for 
the deformed configuration is then 

V, • a = 0, (12) 

in which the Cauchy stress tensor is related to the second Piola-Kirchoff stress tensor in 
Eq. (ll)bycr = /- 1 F-S-F T . 

A finite element simulation using ANSYS 11.0 37 was employed to solve Eq. (12) together 
with the following boundary conditions: 

1. Zero displacements are prescribed at all the side edges of the solid. 
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2. A force balance is prescribed at the top surface. 

n.cr = -p e n, (13) 

where p e is the external pressure. 

3. The pressure at the bottom surface is assumed zero. 

The simulations were carried out for the membrane geometry shown in Fig. 3(a) in the 
absence of a fluid in order to compare the membrane deformation under an external 
pressure loading between a two-dimensional and three-dimensional model. We verified 
with initial simulations of the full geometry mimicking the experimental setup, which 
included the microchannel and pressure chamber, that the deformation was insensitive 
to the physical presence of the pressure chamber and a microchannel, at least for the case 
when the fluid is absent and hence, it was sufficient to simulate a rectangular membrane 
of thickness 60 [im and length 1 mm with fixed edges. The effect of varying the membrane 
width (0.2, 0.5, 1.0, 2.0, 3.0 and 4.0 mm) was examined and compared to a two-dimensional 
finite element model (infinite width assumption (Fig. 3(a))) to determine the limits at 
which the two-dimensional model breaks down. Here, the thin PDMS membrane is 
modelled as a nearly incompressible non-linear neo-Hookean elastic beam with Poisson 
ratio v = 0.495 and Young's modulus E = 2 MPa, as determined from the nanoidentation 
tests in Appendix B. 

III. EXPERIMENTAL DESIGN AND METHODOLOGY 

We fabricated the PDMS microfluidic channel using conventional soft lithography 
processes involving rapid prototyping and replica moulding typically used elsewhere 1 ,2,38 , 
and schematically depicted in Fig. 3(b). To prepare the replica moulds, we spin coat 
multiple layers of SU-8 (SU-8 2035, MicroChem, Newton, MA, USA) negative photoresist 
onto clean silicon wafers, followed by pre-baking on a hot plate at 65°C for 10 min and 
subsequently 95°C for 120 min to remove excess photoresist solvent. Repeated layering 
is required to prepare high aspect ratio moulds with final thicknesses of up to 100-200 
|Um. The photoresist was then exposed to UV radiation at a wavelength of 350-400 nm 
for 60 s through a quartz photomask on which the device designs shown in Fig. 3(c) are 

10 
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FIG. 3. (a) Schematic illustration of the membrane geometry indicating the definition of its width 
and length, (b) Schematic depiction of the main steps involved in the soft lithography procedure 
used to fabricate the PDMS deformable microchannel: (i) deposition of a SU-8 negative photoresist 
layer on a silicon wafer via spin coating, (ii) photoresist exposure to UV radiation through a 
photomask to polymerise the exposed regions, (iii) removal of uncrosslinked photoresist using 
developer solution to produce the replicated structure on the mould, and, (iv) inverse casting of the 
patterned mould in PDMS upon curing, (c) Exploded view schematic showing the microchannel 
(maroon) and pressure chamber (cyan) are cast separately in two PDMS layers and separated by 
an additional thin PDMS layer — the area of this layer which spans the microchannel depth W 
constitutes the flexible membrane. 
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laser printed. This was followed by a two-stage post-exposure bake at 65°C for 1 min and 
95°C for 20 min to enhance the cross-linking in the exposed portions of the SU-8. Finally, 
the wafer was developed to remove the photoresist in developer solution for 20 min and 
the mould dimensions are verified by taking several measurements with a profilometer 

o 

(Veeco Dektak 150, Plainview, NY; 1 A maximum vertical resolution). 

The base and curing agents of two-part PDMS (Sylgard 184, Dow Corning, Midland, 
MI, USA) is mixed in 10:1 ratio and kept in a vacuum chamber to remove any bubbles 
generated during mixing. The PDMS mixture is then poured over the mould and cured in 
an oven at 70°C for 2 hr. To ensure that the rigidity of the PDMS is maintained across the 
devices, the mixing ratio and curing procedure is strictly adhered to. The PDMS channel 
replica was then peeled off the mould and inlet, outlet and pressure sensor ports were 
drilled into the structure. Finally, the PDMS channel was oxidised in a plasma cleaner for 
2 min and sealed by bonding against a flat PDMS substrate. 

The microchannel and pressure chamber are cast in two separate PDMS layers and 
interleaved with an additional thin PDMS layer for the design shown in Fig. 3(c). Mi- 
crochannels with 0.2, 0.5, 1.0, 2.0, 3.0 and 4.0 mm widths were constructed with this tech- 
nique. The thickness of the flexible membrane was fabricated to be 60 fim for all channel 
widths. 

A. Deformation of the PDMS membrane with fluid flow 

The deformation of the thin PDMS membrane under an external air pressure (Preci- 
sion pressure regulator-IR1020, SMC Pneumatics, Australia) applied to the chamber was 
measured visually using a microscope and video camera (AD3713TB Dino-Lite Premier, 
AnMo Electronics, Hainchu, Taiwan; 640 x 480 pixel resolution and 200X maximum mag- 
nification). Nanoindentation experiments, on the other hand, were carried out to evaluate 
the values of the Young's modulus for the PDMS membranes required in the numerical 
simulations, from which values in the range of 1.2 to 2.2 MPa were obtained. A detailed 
description of the nanoindentation test results can be found in Appendix B. Water and 
sucrose syrup were employed as the working fluid, which were driven through the mi- 
crochannel at a constant flow rate using a syringe pump connected to the channel inlet. 
To measure the pressure drop in the channel, we vertically mount capillary tube at the 

12 
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FIG. 4. Images of the patterned SU-8 mould and the PDMS devices for the design depicted in 
Fig. 3(c), panel (i) shows the SU-8 mould for the microchannel, panel (ii) shows that for the pressure 
chamber, and panel (iii) shows a micrograph of the thin interleaving PDMS layer that constitutes 
the flexible membrane that separates the microchannel and the pressure chamber. 



inlet and outlet as illustrated in Fig. 1 and determine the difference in the height across 
the fluid columns in the capillary tubes. 



IV. RESULTS AND DISCUSSION 

A. Comparison Between Two- and Three-Dimensional Models in the Absence of 
Flow 

Figure 5(a) that depicts the maximum deformation predicted by the two-dimensional 
and three-dimensional ANSYS simulations (2D/3D-ANSYS) described in Sec. II B, shows 
that the two-dimensional model begins to deviate from the three-dimensional prediction 
below a membrane width of 2L, when boundary effects associated with edge pinning 
on both sides can no longer be neglected. This is consistent with what we observe in 
the experiments where we measure the deformation of the thin PDMS membrane under 
an external air pressure loading in the absence of fluid flow. Figure 6(a) shows the 
deformed shape of the membrane under an externally applied pressure. The maximum 
deformation of the membrane, measured at the lowest point of the inflexion of its lower 
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FIG. 5. (a) Maximum deformation Az max of the lower surface of the flexible membrane as a function 
of its width W for three different values of the applied external pressure P e/ as predicted by the 
3D-ANSYS simulation described in Sec. II B. The solid lines were added to aid visualisation. Also 
shown by the dashed lines is a 2D-ANSYS (infinite width assumption) finite element simulation, 
from which we observe the diminishing effect of boundary effects associated with the pinned 
lateral edges in the three-dimensional model above a width of approximately 2 mm. (b) Beyond 
this critical width, therefore, we verify that a two-dimensional simulation (2D-FEM-FSI) is a good 
approximation to capture the behaviour of that observed in experiments. 



surface, is extracted visually from similar micrographs and shown in Fig. 6(b) as a function 
of the applied external pressure for microchannels of different widths. The presented 
experimental data point is the statistical average of at least five values, with vertical bars 
indicating the range of the deviation. The experimental error resulted from the manual 
handling of the microscope and video camera, image analysis to extract the deformed 
shape of the membrane and manual measuring of the fluid column height in the capillary 
tubes. In agreement with the predictions of the numerical simulations, we see that the 
deformation becomes independent of the microchannel (and membrane) width when it 
exceeds 2L, therefore suggesting that boundary effects associated with the membrane 
pinning at the lateral edges in a three-dimensional model can be neglected and a two- 
dimensional (infinite width) approximation suffices beyond this critical dimension. 



14 



Fluid-Structure Interaction in Deformable Microchannels 




200 



FIG. 6. (a) Micrographs showing the deformation of the thin flexible PDMS membrane under the 
application of an external pressure p e of 20 kPa for the microchannel designs shown in Fig. 3(b) 
with a channel width of approximately 0.5 mm. (b) Maximum deformation of the membrane 
Az max , measured at the lowest point of the inflexion of its lower surface, as a function of the 
externally applied pressure for microchannels with varying widths. 



The validity of the two-dimensional incompressible neo-Hookean model (i.e., 2D-FEM- 
FSI) is further verified against experimental data for microchannels with large widths 
above 2L. Figure 5(b) shows a comparison between the maximum deformation measured 
in the experiments with that predicted by the two-dimensional model, in that, we observe 
agreement with the experimental data to be bounded by the numerical predictions using 
two values of Young's modulus for the membrane. We note that the large deformation 
data is well predicted by a lower value of the Young's modulus, whereas, better agree- 
ment with the small deformation data is captured using a larger value. Both lower and 
upper values however fall within the 1.2 to 2.2 MPa range measured using the nanoin- 
dentation technique described in Appendix B. In any case, the result provides further 
evidence to suggest that the two-dimensional model is sufficient to capture the membrane 
deformation when it is beyond a critical value of 2L, that is consistently predicted by both 
experiment and simulation. 
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FIG. 7. Relationship between the pressure drop and flow rate across a deformable microchannel 
with two different channel widths for the wider channel. It can be seen that the experimental 
observations (symbols) match well with the predictions afforded by the finite element simulation 
(i.e., 2D-FEM-FSI) described in Sec. II A (dashed line) as well as that for a rigid rectangular 
microchannel (Eq. (15); solid line). For the narrower channel, however, we observe that the 
analytical solution for the rigid microchannel overpredicts the pressure drop at higher flow rates. 



B. Flow Experiments: Pressure Drop and Membrane Deformation 



Figure 7 shows the pressure drop Ap as a function of the flow rate Q obtained from both 
experimental measurements and that predicted by the finite element model (2D-FEM- 
FSI) described in Sec. II A. Also shown is the pressure drop and flow rate relationship 
for two-dimensional fully-developed viscous flow through a long and rigid rectangular 
microchannel, for which, it is possible to obtain an analytical solution if the channel height 
H and width W are small compared to the channel length L c — the solution for longitudinal 
velocity takes the form 39 



Vr = 



4H 2 Ap 
n 3 r]L c 



y - 

^— i n 3 

n,odd 



1 



cosh(nny/H) 
cosh(nnW/2H) 



sin 



tnnz 



(14) 
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Integrating along the width and height of the channel then gives the required pressure 
drop and flow rate relationship: 



_ H 3 WAp 
Q_ Yh\L 



,, n 5 n 



192H , (nnW 



tanh 



.. .z 5 W ^ \ 2H 

i.odd 



(15) 



We observe very good agreement between the experimental pressure drop and flow rate 
relationship and that predicted by both the analytical solution for a rigid microchannel 
and the 2D-FEM-FSI for the case of the wide channel, that lends further support that the 
two-dimensional model constitutes a good approximation when the channel, and hence, 
membrane is sufficiently wide, such that, three-dimensional effects, such as, the pinned 
boundaries at the sidewalls can be neglected. The good agreement with the analytical 
solution, that does not account for the flexible membrane also suggests that the effect 
of the deformation on the pressure drop is small, and hence can be neglected. This is 
however not true for the case for small channel widths where we observe a departure from 
the rigid channel prediction at larger flow rates. Compared to a rigid channel, the ability 
of the membrane to deform in a deformable channel gives rise to smaller effective cross- 
sectional areas that in turn, lead to faster velocities in order that continuity is satisfied. 
Consequently, a lower pressure drop is required to sustain the same flow rate compared 
to a rigid channel with a larger cross-sectional area. 

Figure 8(a) shows profiles of the deformed membrane shape under a specified external 
pressure for varying flow rates (and hence corresponding average inlet velocities (Jo), as 
predicted by the 2D-FEM-FSI simulation described in Sec. II A. We observe no significant 
deformation of the membrane below Uo = 5 x 10" 2 m/s corresponding to a flow rate of 
Q = 110 ml/hr. This is confirmed in our experiments where we do not see any measurable 
changes in the membrane shape at these flow rates. Restrictions in the maximum flow 
rate due to experimental limitations, nevertheless, did not allow us to access flow rate 
regimes where the numerical simulation predicts observable changes in the membrane 
shape. 

Fortunately, however, the membrane deformation is sensitive to the fluid viscosity, as 
shown by the profiles predicted by the numerical simulation in Fig. 8(b). The experi- 
ments were therefore repeated under the same conditions but with sucrose syrup with a 
viscosity of 0.05 Pa • s in order to obtain measurable deformations in the membrane shape. 
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FIG. 8. Shape of the membrane profile as predicted by the finite element model (i.e., 2D-FEM-FSI) 
described in Sec. II A for varying (a) average inlet velocity and (b) fluid viscosity. The rest of 
the parameters are p e = 8 kPa, E - 1.95 MPa, t - 60 jj.m, H - 200 p and W = 3 mm. In (a), 
t\ = 0.001 Pa • s and in (b), U = 5 X TO" 3 . 



Figure 9 shows the close agreement between the shape of the membrane profile that is 
experimentally measured with that predicted by the numerical simulation (i.e., 2D-FEM- 
FSI) described in Sec. II A, therefore inspiring confidence in the predictive capability of 
the proposed model. 

V. CONCLUSIONS 



To investigate the complex fluid-structural interactions between a deformable channel 
wall and the fluid that flows within it, we fabricate a microfluidic device that constitutes a 
single channel out of PDMS with a short section comprising a thin flexible membrane. Ex- 
periments in which we measure the membrane deformation and the pressure drop across 
the channel are complemented by the development of predictive computational models, 
in which, we numerically solve the equations of motion in the fluid coupled with the equi- 
librium equations governing the elasticity of the membrane. In particular, we show that 
two-dimensional models can only be used to describe a three-dimensional system, when 
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FIG. 9. Comparison of the shape of the membrane predicted by the 2D-FEM-FSI (lines) described 
in Sec. II A with that measured experimentally (symbols) for different inlet fluid velocities, (a) 
U = 0, (b) U = 2.2 x 10" 3 , (c) U = 4.4 x 10" 3 m/s and (d) U = 6.6 x 10" 3 m/s. The rest of the 
parameters are p e = 8 kPa, E = 1.95 MPa, t = 60 f/m, H = 200 fim, W = 3 mm and 77 = 0.05 Pa • s. 



the width of the channel, and hence, the membrane is sufficiently large above a critical 
dimension, such that, boundary effects arising from the pinning of the membrane to the 
channel walls at its lateral edges can be neglected — the 2L threshold predicted by the sim- 
ulations agrees well with our experimental observations. In addition, we find excellent 
agreement between the predictions of the deformed membrane shape under an externally 
applied air pressure using both two-dimensional and three-dimensional models with that 
measured in experiments. We believe that the combination of these results, the predic- 
tive capability of the numerical models developed, and the physical insight gleaned in 
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this study would be useful in the design of polymer-based microfluidic devices, and, 
in particular, microactuation schemes such as the pneumatically-driven micropumps, 
micromixers, microvalves and microfilters employing flexible polymer membranes that 
have grown increasingly popular over the last decade. 
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Appendix A: Weighted Residual Form of the Equilibrium Equation Vx • S = 

Here, we provide a correction to the weighted residual form of the equilibrium equation 
given by Eq. (4) derived by Carvalho and Scriven 29 . The error in the original derivation, 
whilst insignificant for small deformations, leads to significant discrepancies when the 
deformation is large. 

The weak form of Eq. (4) is 

f (V x • S)<(> dn s = - f (V x • S) dQ s + f <MN • S) df s = 0, (Al) 

Jq s Jq s Jt s 

where Qs, Ts and N are the area, arc length and unit normal in the zero-stress configu- 
ration, respectively, and (p is a weighting function. When written in terms of Cartesian 
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components, the weighted residual form of this equation in the computational domain is 
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Here, Qso and T$o is the area and arc length in the computational domain, respectively, |J*| is 
the Jacobian of the transformation from the zero-stress configuration to the computational 
domain and (pi is a bi-quadratic weighting function. The components of the dimensional 
Piola-Kirchhoff stress tensor S in terms of the dimensional Cauchy stress tensor for a 
neo-Hookean material a = -n* I + G B, are 



c JV ^ r dx Jy dx 
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where I is the identity tensor, n* is a pressure-like scalar function and B is the left Cauchy- 
Green tensor. In their finite-element formulation of the fluid-structure interaction prob- 
lem, Carvalho and Scriven 29 (see also Carvalho 40 ) have used 
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(A6) 



in place of Eqs. (A2) and (A3). Essentially, the positions of the two components Sy x and 
Sxy have been interchanged. 

In order to establish the validity of Eqs. (A2) and (A3) and to demonstrate the incorrect- 
ness of equations (A5) and (A6), we have examined the simple problem of a beam fixed 
at the edges with uniform pressure applied on both the top and bottom of the beam, as 
shown schematically in Fig. 10. Essentially, we compare the results of our computations 
using Eqs. (A2) and (A3) (labelled FEM-N), and Eqs. (A5) and (A6) (labelled FEM-C), with 
the results obtained with the finite element ANSYS simulation for a plain-strain model 
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FIG. 10. Geometry of the solid domain. 

described Sec. II B. In addition, we prescribe boundary conditions in the form of zero 
displacements at the left and right edges of the beam, and a force balance at the top and 
bottom of the form, 

n-cr = -p i n, (1 = 1,2) (A7) 

where n is the unit normal to the deformed solid surface, and p\ and p2 are the dimensional 
external pressures on the top and bottom of the beam, respectively. 

In units of height H, the length of the beam is set at L = 5H, with H = 10~ 3 m. The 
external pressures have been chosen to be pi = 1.1 N and p2 = 1.0 N, and three different 
values (6000, 12000 and 24000 Pa) have been used for the shear modulus G. Computations 
have been performed with three different meshes (Ml, M2 and M3) in order to examine 
mesh convergence. We note that the formulation of the fluid-structure interaction problem 
by Carvalho and Scriven 29 applies to the special case of an incompressible neo-Hookean 
material with a Poisson ratio v = 0.5. On the other hand, the ANSYS plain-strain package 
is only applicable to compressible neo-Hookean materials. Consequently, in order to carry 
out the comparison with the ANSYS simulations, we have obtained predictions with 
several values of v < 0.5, and extrapolated the results to v = 0.5. 

Figure 11 compares the maximum displacement of the beam obtained with the FEM-N 
and FEM-C formulations with the ANSYS plain strain model for the three different values 
of G. In all three cases, mesh-converged results obtained with FEM-N are observed to 
agree with the extrapolated mesh converged solution obtained with ANSYS. On the 
other hand, the mesh converged solution obtained with FEM-C shows differences from 

22 



Fluid-Structure Interaction in Deformable Microchannels 



2.10e-004 



* 




& 1.80e-004 



1.50e-004 



G=24000 Pa 







0.005 



O 
* 


+ 

v 
x 
* 
D 

0.01 0.015 
0.5-v 

(a) 



ANSYS Ml 
ANSYS M2 
ANSYS M3 
FEM-NM1 
FEM-N M2 
FEM-N M3 
FEM-C Ml 
FEM-C M2 
FEM-C M3 



0.02 



0.025 



4.50e-004 



& 3.25e-004 

■■Z 



G=12000 Pa 







2.00e-004 



o 



+ 

v 
x 

* 
D 



ANSYS Ml 
ANSYS M2 
ANSYS M3 
FEM-N Ml 
FEM-N M2 
FEM-N M3 
FEM-C Ml 
FEM-C M2 
FEM-C M3 







0.005 



0.01 



0.015 
0.5-v 

(b) 



0.02 



0.025 



6.55e-004 



3 



1 6.25e-004* - 



3 



g 5.95e-004 



G=6000 Pa 



5.65e-004 




■ ANSYS Ml 
O ANSYS M2 
* ANSYS M3 
FEM-N Ml 
+ FEM-N M2 
V FEM-N M3 



0.005 



0.01 0.015 
0.5-v 



0.02 



0.025 



(C) 

FIG. 11. Comparison of formulations FEM-N and FEM-C with the ANSYS plain-strain model for 
three different values of G: (a) 24000 Pa, (b) 12000 Pa, and (c) 6000 Pa. 
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the other two approaches. Indeed, while this difference is small at G = 24000 Pa, and 
substantially larger at G = 12000 Pa, we are unable to obtain a converged solution with 
FEM-C at G = 6000 Pa. 

Appendix B: Nanoindentation Characterisation of the Elastic Modulus of the PDMS 
membrane 

Since the Young's modulus of PDMS can significantly be altered by varying the cur- 
ing temperature and time, as well as the mixing ratio of the silicone base to the curing 
agent 2,8 " 10 , there is a need to measure the isotropic mechanical properties of PDMS 41 " 43 . 
Different experimental techniques have been employed to characterise the rigidity of 
PDMS and the reported value of Young's modulus for PDMS usually falls within the 
range of 0.05-4.0 MPa 8,9 . Recently, Liu et alA 1 have conducted a tensile test to estab- 
lish the thickness-dependent hardness and the Young's modulus of PDMS membranes, 
arising due to the shear stresses that are exerted during fabrication of these thin mem- 
branes. On the other hand, nanoindentation testing, which has been widely employed 
for characterising the elastic and plastic properties of hard materials, is also now being 
recognised as a tool for characterising the mechanical properties of polymeric materials. 
In a standard nanoindentation test, the Young's modulus and hardness of a very thin 
membrane made of elastic material can easily be obtained from the load displacement 
data. Carrillo et a/. 44 , for example, used a nanoindentation technique to characterise the 
Young's modulus of PDMS with different degrees of crosslinking. 

Here, PDMS (Dow and Corning Sylgard 184) samples were prepared by mixing the 
cross-linker and siloxane in a ratio of 1:10, and subsequently kept in a vacuum chamber to 
remove the bubbles that were generated during mixing. PDMS membranes with different 
thicknesses were then produced by spin coating glass wafers at various rotation speeds 
followed by curing in an oven at 70°C for 2 hrs. The thickness of the PDMS membrane 
was measured using a surface profiler. By varying the rotation speed of the spin coating 
process between 500 and 2000 rpm, PDMS membranes of thicknesses in the range of 25[im 
to 100 /im were produced. 

The nanoindentation testing is carried out using a Tribolndenter® (Hysitron, Inc., Min- 
neapolis, MN, USA) at room temperature with a Berkovich indenter tip. For load control 
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function, we employ a loading and unloading rate of 10 /iN/s, peak load of 100 (jN and a 
hold period of 5 s. When the tip of the indentor reaches the sample surface, the instrument 
applies the predefined load and records the load and displacement data accordingly. The 
hardness and the Young's modulus of the material is then determined from the unloading 
portion of the load-displacement curve using classical Hertzian contact theory 45 : 

H = %^, (Bl) 

A 

K = — + ^P <B2) 

where H is the hardness of the substrate and F max is the maximum force applied on the 
PDMS membrane. A is the projected contact area between the tip and the substrate, v 
and E are the Poisson's ratio and the Young's modulus for the test specimen, respectively, 
and V; and E, are those for the indenter. The material properties of the diamond indenter 
are E, = 1140 GPa and v,- = 0.07. The reduced elastic modulus E r is calculated using the 
following expression proposed by Oliver and Pharr 46 : 

- SV5 (B3) 



Ifi^A 

where S is the contact stiffness, taken as the initial slope of the unloading section of the 
load-displacement curve, and /3 is a constant that depends on the geometry of the indenter. 

Figure 12(a) shows the indentor displacement in response to the applied load during 
a nanoindentation test carried out on the PDMS membrane using quasi-static measure- 
ments. To confirm the reproducibility of the test data, the indentation is performed on 
nine different locations of a given sample. The penetration depth of the indenter is ob- 
served to be higher for low thickness PDMS membranes indicating a lower value for 
Young's modulus as the thickness decreases. Similar trends are also observed for other 
PDMS membranes. 

The reduced elastic modulus E r and Young's modulus of the PDMS membrane are cal- 
culated using Eqs. (B2) and (B3). Figure 12(b) indicates that the Young's modulus of PDMS 
membranes decrease with decreasing membrane thicknesses. This is due to the polymer 

25 



Fluid-Structure Interaction in Deformable Microchannels 




500 1000 1500 2000 2500 3000 3500 4000 4500 
Displacement(nm) 



3.5 

£2,5- 
1 

a i.5- 
■8 

a 

a if 

0.5 





T • 

(I 

A 

J I 



0.05 
Thickness(mm) 



0.1 



(a) 



(b) 



FIG. 12. (a) Load-displacement curves for a PDMS membrane with a thickness of 0.05 mm. (b) 
Thickness dependence of the reduced elastic modulus E r and Young's modulus E of the PDMS 
membrane. 



molecules experiencing enhanced radial stretching due to an increase in the shear stress 
with the increasing spinning speeds that are required to produce thinner membranes. In 
any case, the good agreement between the measured values for the Young's modulus 
therefore suggest that the nanoindentation test is capable of differentiating the elastic 
behaviour of a polymeric material with varying thickness. 

Appendix C: Validation of the Finite-Element Formulation 

Here, we briefly provide results on the validation of the finite element formulation 
described in Sec. II B against simple benchmark cases that have been reported earlier in 
the literature. 

1. Couette Flow Past a Finite Thickness Solid 

The flow of a Newtonian fluid past an incompressible neo-Hookean solid, as shown 
schematically in Fig. 13, has been previously described by Gkanis and Kumar 47 . The 
interface between the fluid and solid is located at y = t and a rigid plate located at z 
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FIG. 13. Schematic depiction of Couette flow of a Newtonian fluid past an incompressible finite 
thickness neo-Hookean solid. 



= (H + t) moves in the x direction at a constant speed U , giving rise to Couette flow 
in the fluid domain; the bottom edge of the solid is held fixed. Gkanis and Kumar 47 
performed a linear stability analysis of this problem in the limit of zero Reynolds number 
Re and infinite domain length L, and have shown that the steady-state solution of the 
deformation in the solid produced by the Couette flow is 



x = (x,y) = (X + TY,Y), 



(CI) 



where T is the dimensionless number defined as T = (i]Uo)/(GH). 



Computations have been performed to compare predictions for the deformation of the 
solid domain at L/2 with the analytical results of Gkanis and Kumar 47 . In order to elimi- 
nate end effects caused by the fixed ends of the solid and fluid domains in computations, 
we have varied the length of the domain between 10 and 30 m and ensured that domain 
length independent predictions are obtained. The following parameter values have been 
used: p = 10" 3 kg/m 3 , /] = lPa-s,H=lm,f = lm and U = 1 x 10" 3 -1. 75 x 10" 3 (such 
that Re ~ 0) and G = 10" 2 Pa. This choice of parameter values maintains T in the range 
0.1-0.175. The mid-surface displacement of the solid predicted by the finite-element for- 
mulation is compared with the analytical solution for different values of T. Figure 14 
shows that in all cases the predictions of our finite element simulation are in excellent 
agreement with the analytical solution. 
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values of T. The lines denote the analytical solution reported by Gkanis and Kumar 47 whereas 
symbols are the predictions from the our simulation. 



2. Flow in Two-Dimensional Collapsible Channels: Elastic Beam Model 

Luo et al. 48 have carried out extensive studies of Newtonian fluid flow in a two- 
dimensional collapsible channel by considering the flexible wall to be a plane-strained 
elastic beam that obeys Hooke's law. In contrast to the current finite thickness elastic solid 
model, the beam model does not admit any stress variation across the beam cross-section. 

For the purposes of comparison, the dimensions of the channel and other parameter 
values are chosen to be identical to those used by Luo et ah 48 in their simulations: L u = 5H, 
L = 5H, U = 30H, U = 0.03 m/s, H = 10" 2 m, p = 10 3 kg/m 3 and r\ = 0.001 Pa • s. This 
choice corresponds to Re = 300. Further, we set G = 11.97 kPa (which is equivalent to 
a value of 35.9 kPa for the Young's modulus of an incompressible solid) and p e = 1.755 
Pa. The flexible wall thickness is varied in the range 0.01H-0.1H. We note that the 'pre- 
tension' in the beam is also a variable in the model of Luo et al. 48 ; however, since no 
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FIG. 15. Comparison of the shape of the flexible wall predicted by the finite thickness elastic solid 
model with the results of Luo et al. 48 for an elastic beam model with different values of the wall 
thickness. Circles denote the elastic beam model whereas lines denote the results of our numerical 
simulation. 



such variable exists in the current model, we have restricted our comparison to the results 
reported by Luo et al. 48 for cases where the bream pre-tension is zero. 

Figure 15 compares the prediction of the shape of the flexible wall of our finite thickness 
elastic solid model with that reported by of Luo et a/. 48 . While our simulations agree 
with Luo et al. 48 for the relatively small deformations that occur at large membrane 
thicknesses t, the Hookean beam model, as expected, begins to depart from the prediction 
of the nonlinear neo-Hookean model for large deformations that are associated with small 
membrane thicknesses. 
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3. Flow in Two-Dimensional Collapsible Channels: Zero-Thickness Membrane 
Model 

Simulations have also been performed to compare predictions of the flexible wall shape 
by current finite thickness elastic solid model with that of a zero-thickness membrane 
model of Luo and Pedley 23 for the flow of a Newtonian fluid. Apart from the simplicity 
of the zero-thickness membrane model from a constitutive point of view, a fundamental 
difference between the two models is that while the tension in the flexible wall is prescribed 
a priori in the zero-thickness membrane model, it is part of the solution in the finite 
thickness elastic solid model. As a result, a several step procedure is required to carry out 
the comparison, as described in what is to follow. 

The zero-thickness membrane model is first computed for a pre-determined value 
of membrane tension equal to 675 N/m, with the following parameter values: Re = 1, 
p = 1054 kg/m 3 , U = 1.338 x 10" 2 m/s, H = 10" 2 m, r\ = 0.141 Pa • s and p e = 17545 N/m 2 . 
This leads to a prediction of the minimum height of the gap in the channel (beneath the 
flexible membrane) of h/H = 0.125. Computations with the finite thickness elastic solid 
model are then carried out for the same parameter values, for various combinations of 
flexible wall thickness t and shear modulus G such that each combination always leads 
to the same value of the minimum channel gap height, namely h/H = 0.125. It turns 
out that even though the minimum gap height is the same in both models, the predicted 
interface shape is not, with the difference increasing as the thickness of the elastic solid 
increases. This is clearly a result of the finite thickness of the elastic solid. Consequently, 
in order to compare the interface shape, we carry out an extrapolation procedure in which 
the interfacial height of the interface at various locations in the gap as a function of the 
flexible wall thickness is extrapolated to the limit of zero wall thickness, as shown in 
Fig. 16(a). The extrapolated interface shape is then compared with the prediction by the 
zero-thickness membrane model in Fig. 16(b), in which we observe excellent agreement 
between the two models. 

The remaining step involves the evaluation of the resultant tension in the finite thick- 
ness elastic solid and how it compares with the pre-determined membrane tension of 
675 N/m. Here, we first estimate the tension in the finite thickness solid at a particular 
location x by averaging the tangential solid stresses acting across the cross-section at x. 
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FIG. 16. (a) Extrapolation to t — of the flexible wall shape obtained from the finite thickness 
elastic solid model for t — 0.01H, t — 0.05H, and t = 0.1H. (b) Comparison of the shape of the 
flexible wall predicted by the finite-thickness solid model (symbols) with the prediction of the 
zero-thickness membrane model (solid line), (c) Extrapolation of the average tension acting in the 



elastic solid to the limit of zero wall thickness. 
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An estimate of the overall tension in the solid is then obtained by averaging the tension 
along the entire length of the flexible solid for all values of x. The values of the average 
tension obtained from the finite thickness elastic solid model for t = 0.01H, t = 0.05H 
and t = 0.1H are then extrapolated to t = 0, as shown in Fig. 16(c), in which we observe 
the extrapolated value of tension to be fairly close to the value of 675 N/m used in the 
zero-thickness membrane model. 
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